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Abstract 

The studies of Dyson-Schwinger Equations (DSEs) provide us with insights into nonperturbative phenomenon 
of quantum field theory. However, DSEs are essentially an infinite set of coupled Green’s functions, it’s necessary 
to decouple parts of the equations which are thought of major physical importance to make the solution of these 
equations possible. Although the results are model-dependent, no qualitative deviations from exact solutions are 
expected with properly chosen truncation scheme. In this article, a globally convergent numerical method for 
the solution of the DSEs of QED 3 in Euclidean space is presented. This method can be adapted for more complex 
problems, however, it also shows its limitations when adopted in problems such as the searching for Wigner solutions. 

1 Introduction 

Quantum Electrodynamics in ( 2 + 1 ) dimensions (QED 3 ) is studied for a variety of reasons pTflQ]. It possesses 
nonperturbative features such as confinement and dynamical mass generation while avoids additional complication 
of non-Abelian theory HU. Eurthermore, it has also been proposed that QED 3 can be employed to study the High-Tc 
cuprate superconductors, where the motion of electrons are confined to two-dimensional planes HMl]. However, 
as the analysis of nonperturbative solution cannot be conducted on an order-by-order basis from certain starting 
point, the solution of DSEs relies on self-consistent physically motivated models HU. 

Although the normal fixed-point iteration works well in the calculation DSEs with bare vertex and the minimal 
Ball-Chin (BCl) vertex [16], this method can’t provide sufficient sampling at the the singular point when calculating 
the dressing function using the full Ball-Chin (BC) vertex [T7]. This problem can be resolved by adopting the 
polynomial expansion of of the unknown function [18] and apply Newton Method to the iteration of the coefficients. 
In this method, as long as coefficients of the expansion are known, we can calculate the unknown function over 
the range of expansion. The infrared and ultraviolet behavior can be determined analytically. However, the local 
Newton depend critically on the initial values of iteration, therefore a algorithm with global convergence is adopted 
to reduce the dependence on the initial value of iterations. 


2 Formulation of Dyson-Schwinger equations in QED3 

2.1 functional formulation of QED3 

The Lagrangian density of QED 3 with Nf fermion flavors can be written as 


= -ml +el4)ipf - 

/=i ^ 

where / is the flavor index, Nf is the total number of flavors and ^ is the gauge-fixing constant. 

The generating functional of correlation functions defined via this Lagrangian density reads: 

Z[ri,fj,Ji_,] = J[di)4x)][di^Hx)]d[A^(x)]exp\i jd^x[^{i},%l),A^) + ^{ii)fr]^+rif%l)f) + A^J^] 


The generating functional of connected Green’s functions rj, J^] is related to Z[rj^ rj, J^j] by HS] 

z[n, fj, Jf,] = exp{Gc[n, fj, J^]). 


( 1 ) 


( 2 ) 

( 3 ) 
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We can perform the Legendre transformation of Gc[r]^ J^] to obtain the generating functional of proper vertices 


iT[ip,'ip,A^] = - i J (f x[tp^7]^ + rj^ 


with identifications 

/ f _ 7f_ 

~ iSfjf^ ~~i5i^r 


A,= 


SGe 

iSJ^ 


( 4 ) 


( 5 ) 


It can be shown that r[^, ip, Afj\ generates the 1-PI Green’s functions [Hj. Then the fermion propagator S{x, y), 
photon propagator and and the 1-PI fermion-photon vertex can be expressed as 


iS^{x,y) = — 


S^Gc 


iDf,„{x,y) = - 


iSrjf (x) (y) 

S^Gr, 


77=0,77=0, J =0 


elT{x,y;zy = 


i5J>^(x) i5J''{y) 
5^T 


77=0,77=0, J =0 


5Ayz)57p{x)5%p{y) 


77=0,77=0, J =0 


For simplicity, the flavor index will be omitted in the following context. 

2.2 Dyson-Schwinger equations for the fermion and photon propagators 

The complete two-point Green’s function S can be constructed out of the electron self-energy Sp by 

iS{p) = iSoip) + iSo{p)[-i^{p)]iSo{p) + i5ob)[-i5](p)]i5ob)[-i5](p)]i5ob) + • • • 

where 

1 


Soip) = 

then we have 
Sip) = 


— rtiQ ’ 

1 


^ - Too - E(p)' 


(6) 

( 7 ) 

(8) 


( 9 ) 

( 10 ) 

( 11 ) 



*> 
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Figure 1: The complete fermion propagator as a sum of IPI self-energy insertions. 

The Schwinger-Dyson equation for the fermion self-energy is given by 

/ (Pk 

j^^lniS{k)Vyk,p)iD^yq), ( 12 ) 

where q = p — k. 

Similarly, the complete photon propagator can be written in terms of self-energy as 

W^yq) = W^yq) + iD^yq)[iU,x{q)]iD^''{q) + ■ ■ ■ (13) 

where 

^o'^(9) = ^ + (i-e^, (14) 

we obtain 

D^yq) = D^yq) - D^yq)U,x{q)D^yq) + ■ ■ ■ (15) 
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Figure 2: Diagrammatic representation of the integral equation for fermion self-energy S(p). 



+ 



Figure 3: The complete photon propagator as a sum of proper IPI photon self-energy. 


The DSEs for the photon self-energy reads 

/ d^h 

^Tr[iS{k)^^iS{k + q)V,{k + g, k)]. (16) 



Figure 4: Diagrammatic representation of the integral equation for 


2.3 Dyson-Schwinger equation for the fermion-photon vertex 

The Dyson-Schwinger equation for fermion-photon vertex can be shown to be m 

/ d^k 

+ k)Tf,{p' + k,p + k))iS{p + k)K{p + k,p' + fc, k), (17) 

where K is the fermion-antifermion scattering kernel. 

The field theory can be completely determined when all of its Green’s functions are known. DSEs connected 
all these n—point functions through a set of infinitely coupled integral function. In the derivation for the fermion- 
photon vertex, we encounter another unknown 4-point function K which satisfies its own integral equation and 
will be coupled to a 5-point Green’s function.To make the DSEs tractable, it necessary to introduce the truncation 
at certain n. We can choose as the building block of our theory. Quantities such as K be expressed 

in terms of by the analysis of their topological structure. The theory can then be determined self- 

consistently in this way. To compensate for the strategy of truncation, we have to introduce physically motivated 
vertex in our calculation. Common restriction on vertex are various symmetries and the known asymptotic behaviors 
of dressing functions. Several attempts have been made to determine the form of Ty [22j[24]. The simplest vertex 






is the bare one, however, this violate the Ward-Takahashi identity which is a consequence of the gauge invariance 
of the theory. In this article, a numerical method for the solution of DSEs in the BC vertex which satisfied the 
Ward-Takahahsi identity is presented. Due to the singular structure of the kernel of integration, the fixed point 
can’t be to calculate the dressing function with precision to satisfactory degree. The Newton’s iterative method is 
adopted in the literature [18]. However, the dressing functions are not calculated in the same loop in the method 
proposed by the author [18] . And the author has to manually cancel the infrared divergence which is not supposed to 
appear in the final result of calculation. The method presented in this article will solve the three coupled equations 
consistently. 








Figure 5: Diagrammatic notation for the Dyson-Schwinger equation of fermion-photon vertex. 


2.4 The truncation scheme 

Following the construction in [24], it can be shown that the BC vertex takes the following form 
A(p^^A(k^)) (A(p^)-A(k^)), B(p^)-B(k^), 


2(p2 — k‘^) 


— k‘^ 


(18) 


To perform the integration in 3-dimensional Euclidean space, we transcribe the Dyson-Schwinger equation by 
the following rules 


o 

1 

t 

o 

(19) 

X —>• X, 

(20) 

^ -ik°, 

(21) 

k^k, 

(22) 


A 





















and 


7f =7^ 


-27 


(23) 

(24) 


So{p) = - 

Due to the spinor structure of fermion propagator, it admits the following general form 
S~^{p) = ii)A{p^) + B{p^). 

Substitute Eq. (EZI into DSEs, we obtain the following set of couple integral equations 

(fk 




li 

so that 

{7^,74 = 2V- (25) 

Note that we also need to perform a reflection of k so that the free-fermion propagator has the form 

^ (26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 


B{P^) = \ I 


^TrMk)r.D,.{q)]. 


In Euclidean space, the photon propagator has the form 
d^k 


n„ 


A unit in which = 1 has been chosen, it can be written 

n{qY^ = - q^n^iq^), 

by the Ward-Takahashi identity 

g^n^"(</) = 0 . 

The ultraviolet divergence in n(g^) can be projected out by the following operator [ 2 T] 
qv 


.Tr[S{k)-f^S{p)T,{p,k)]^ 


Substitute the vertex into Eqs. (EH), (EH) and (Eni and employ the trace techniques list in Appendix A, we arrive 


at 


A(p 2 ) =1 


+ 
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d?k f IA(pY + A(kY]A(kYl(p^k^ - ip ■ kY)/q^ -p-k] 


(27r)3 V 


[A2()!c2)fc2 B2()!c2)]g2(l n(g 2 )) 


[2(p2fc2 _ (p . kf)]{B{p^) - B{k^))B{k^) 

{A^{k^) + B2(fc2))^4(l + n( 92 ))(p 2 _ k2) 

A{kY[AipY - A{kY] 


b{pY = I 


+ [ 2 (A:^ +p 2 )(A:V - ip ■ kf)] 

d^k 1 


(27r)3 {A^{k^)k^ + B2(fc2))(i + n(g2)) 
X [ 2 { p ^ k ^-{ p - kY )/ q ^)]}, 


2{A^{k^)k^ + B^{k^)){p^ - fc2)g4(i 4 . n(g 2 )) / ’ 

f BjkYiAjpY + A{P)) B{P)AipY - AjkYBjpY 


n(5^) 

+ 


d^k 


A{pYA{kY[A{p‘^) + A{kY] 


(27r)^ {A^{k‘^)k‘^ + B‘^{k‘^)){A‘^{p‘^)p‘^ + B‘^{p‘^)) 
{A{p^) - A{k^)) 


q 2 (^ p 2 _ ^2^ 

[2k^ -4:{k-q) -6{k-qY/q'^] 


(34) 


(35) 


[4(B(fc2)B(p2))((^2 ^ 4 ^ . ^ ^ . qY/qZ _ 


2{A^(kP‘)kP‘ + B(}P‘)){A^{p‘^)p‘^ + B{p‘^))(jA — kP‘) 

+ 2A{k^)A{p^){Ak'^ - Ak^k • q - 10{k • qf - 2k • qq‘^ - 12k‘^{k • qf/q‘^ - 12{k • qf/q‘^)] 
+ 4:[B{k‘^)A{p‘^){k‘^ — 3{k • q)^lq‘^ — 3k • q — q‘^) B{p^)A{k^){k^ — k ’ q — 3{k • g)^/g^)] 

B{p‘^) - B{k‘^) 

— kP‘){A^ {k‘^)k‘^ + B{k‘^)){A^{p‘^)p‘^ + B{p‘^)) 


}■ 


( 36 ) 


Our task is to find the numerical solution of these equations. 


















3 Numerical method 


3.1 Polynomial expansion 

As an improvement over the fixed-point integration, the polynomial expansion of the unknown function is 
adopted. This strategy is necessary if we want to calculate the values of function at arbitrary points (The fixed-point 
integration, as indicated by its name, only calculate values of functions at fixed point, and reduce the continuous 
values of functions to a discrete set of function values. In the calculation of complex vertex such as BC vertex, we will 
have to handle the singular integral kernel. We would encounter singular term proportional to l/{p‘^ — We may 
feel attempted to skip over the singular point of integration, however, this would cause the loss of accuracy and the 
iteration may never converge. On the other hand, if we expand the unknown function in terms of known polynomials- 
for example, the Chebychef polynomial-as long as the coefficients of expansion is known, we can calculate the value 
of the functions at any points. In this methods the coefficients will be calculate iteratively in Newton’s method). 
Chebychef polynomials have the following form 


Ti{x) = cos{i acos{x)), 

(37) 

where x G [—1,1]. 

We also need to change the variable of integration to ln{p^) where p is the momentum, to 
sampling at singular points. Numerical cutoff is also needed so that the integration is done over 
Therefore we arrive at the mapping 

ensure sufficient 
the range [ 5 , A]. 

ln(P) 

2 

(38) 

In general the unknown function can be expanded as 


Mp^) = Y + 

(39) 

B{p^) = ^ + El,biTiix), 

(40) 

n(p2) = 1 + E^.aTiix), 

(41) 


where N is the number of terms in the approximation and Ti is the i th term of the Chebychef polynomial. 


3.2 Newton’s method 

The Newton iteration 

Xn+l =Xn- F'{Xn)~^F{Xn), 


(42) 


exhibits q-quadratical convergence [23]. Therefore it was chosen as the basis of the algorithm. However, we would 
have to deal with the addition complication of the derivatives of the unknown functions with respect to the coefficients 
of the expansions. The integration in Eqs. (|34|) . (|35)) and (|36|) can expressed in the spherical coordinate in 3- 
dimensional Euclidean space 

2 \ 1 /* 1 


• ^.^([A{p‘^)-\-A{k‘^)]A{k^)[{p‘^k‘^sin‘^0)/q‘^—pkcosO] 
sinOdO s - tfz - ^' - 


B{P^) = I 


(p2-fc2)(i + n(g2)) 

k'^dk 1 


2 (p 2 -fc 2 )(i+n(g 2 ))' 


(27r)2 (242(fc2)fc2 +52(^2)) 


I s. 


sinOdO 


1 


(i + n(g 2 )\- 


r B{k^){A{p^) + .4(fc2)) B{k^)A{p^) - A{P)B{p^) 


[ 2 p 2 fc 2 




k^dk 


(27r)2 [A^{k‘^)k‘^ + B{k'^)) 




sin OdO 


X [2k‘^ — AkqcosO — 6{kq cosO)^/ 


W) - A{k^)) 


— k‘^) 

{A{p^)A{e)[A{p^) + A{e)] 


}' 


(43) 


(44) 


— k‘^) 


[2{B(k^)B{p ))(g + AkqcosO -|- 6{kq cosOY/q 


+ A{k‘^)A{p‘^){4:k^ — Ak^kqcosO — l{){kqcosOY — I2k^ {kq cosOY / q^ — Vl(kqcosdY j Y ~ 2(kqcosS)Yy\ 


+ 


B{p^) - B{e) 


4:[B{k‘^)A{p‘^){k‘^ — 3{kq cosOY / q‘^ — 3kqcos0 — q^) 


(p^ — k^) 

+ B{p^)A{k^){k^ — kqcosO — 3{k qcosOY/q‘^)]} , 




2Y) 


( 45 ) 

















Substitute the expansions (l39l) , (I41II and (HOll , we obtain the following nonlinear system of integral equation 


FA{ai,bi,Ci) = 0 , 
FB{ai,bi,Ci) = 0 , 
Fn{cii,bi, Ci) = 0 , 


(46) 

(47) 

(48) 


where 


FA{ai, bi, a) =A{p^) - 1 + 


?/ 


k^dk 


■Jsi 


sinOdO 


[A{p^) ^ A{k^)]A{e)[{p^e sin^0)lq^ -pkcosO] 


p‘^ J ( 27 r)^ {A^{k‘^)k?‘ + B‘^{k‘^)) 

{B{p^) - B{k^))B{e) , ^(fc2)[^(p2)_^(fc2)l ^ 


52(l + n(g 2 )) 


' ' (p2 - fc2)(l + n(g2)) ' 2(p2 _ + n(g2)) + ^^ )] [2(^ p ^)/<l 

/ k^dk 1 

(27r)2 {A^{k^)k^ + B^k^))^ 

1 r B(fc2)(^(p2) + A(fc2)) B{P)A{p^) - A{k^)B{p^) 


(49) 


/ sinOdO- -< 

(i + n(</ 2 ) \ 

r . ^ TT/ 2x Nf f Pdk 1 

n(ai, i,Ci)- (p ) ^2 y (27r)2 (^ 2 (fc 2 ))^ 2 + 5 ()^ 2 )) 


/ 


+ 


sin OdO 


q^{p‘^ — kP‘) 

1 

A^{p‘^)p‘^ ^ B{p^) 


'}■ 


[ 2 p 2 fc 2 

^ (50) 

{A{p^)A{e)[A{p^) + Aie)] 


2 , 0 , , {A{p‘^)-A{e)) 


X \^k^ — 4:kq cosO — ^ikq cosOY/q ] ^ 

[p^ — k^) 

X \^{B(jY)B{p^)){q^ + AkqcosO + 6(/cg' cosS)^j— 2k‘^) 

+ — Ak^kqcosO — l{){kqcosO)‘^ — 12k^{kqcos0)‘^ / q^ — \2(kqcosdkf j — 2(kqcosG)c^y\ 


^^ 4:[B{k‘^)A{p‘^){k‘^ — ‘^{kqcosO)‘^/q^ — SkqcosO — q^) 
[p^ — k^) 

B{p^)A{k^){k^ — kqcosO — 3{k q cosO)‘^/q‘^)]} , 


(51) 


and 

A{p^) = ^+El,aiTi{x), 

B{p^) = ^-^ + El,biTiix), 
n(p 2 ) = I + s(I,CiTi(x). 

We next calculate the derivatives of FA{cii^bi^Ci)^ Fsiai^bi^Ci) and Fuidi^bi^Ci) with respect to ai^bi and q, 


respectively. To clarify the calculation, the following notations are adopted 

Ti = {p‘^k‘^ sirY0)/q‘^ —pkcosO, (52) 

T 2 = 2{k‘^p‘^ sirY0)/q^, (53) 

Ts = 2k^ — AkqcosO — 6{kq cosO)^ / q^, (54) 

T 4 = + AkqcosO + 6{kq cosG)^/cY — 2k^^ (55) 

T 5 = FdY — Ak^kqcosO — l^{kqcosO)^ — I2k^ {kq cosO)^ / q^ — 12 (kq cosO)^ / q^ — 2{kqcos0)q^, (56) 

Tq = k^ — 3{kq cosO)^ / q^ — 3kqcos0 — g^, (57) 

Tj = kY — kqcosO — 3{kqcos0)^/q^ ^ (58) 
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then we have 


OFa 

dtti 


=Ti{x 


+ 


>-/ 


k‘^sinO dk dO 


A{e) 


(27r)^ + Ii{q^)){A^ {k‘^)k‘^ + B‘^{k‘^)) 

A{k‘^){jp‘ + 


TiTi{x) 


2{p^ — k‘^)p^q^{l + Ii{q^)){A^ {k‘^)k‘^ + B‘^{k^)) 

2A^{e){A{e) ^ A{p^))e 


T2Ti{x) 


TiTiiy) - 


+ 


p 2 g 2 (i ^ n(g 2 ))(A 2 (A: 2 )p + B‘^{k‘^)f 

A^{k‘^){A{p‘^) — A(k^)){p^ + k‘^)k^ 

{p‘^ — k‘^)p‘^q^{l + Il{q^)){A^ {k‘^)k‘^ + B‘^{k‘^)Y 
{A{P) + A{p^)) _ 

p 2 g 2 (l + n(g 2 ))(^ 2 (fc 2 )p + 52 (fc 2 )) 1 AU) 

A{p^) — A{k'^) 

2 (p 2 — k?)p^q^{\ + 'n.{q^)){A'^ {k‘^)k‘^ + B^{k^)) 


2A(fc2)B(fc2)(E(p2) - B{k^))P 
(p 2 — k^)p^q‘^{l + Il{q'^)){A‘^{k^)k^ + B^{k^))^ 
A{k‘^) 


T2Ti{y) 


T2Ti{y) + 


TiT^iy) 

T2Ti{y) 


p^q^il + n(g 2 ))(A 2 (A: 2 )fc 2 + ^ 2 (^ 2 )) 

_ A{P){p‘^ + e) _ 

2(p2 — k'^)p'^q^{l + n((/2))(^2^^2j^2 _|_ q2Q,2^'^ 

T 2 Ti(y)|, (59) 


OFa 

dbi 


f k'^sinOdkdO f 

J (27r)2 \ 


B{P) 


(27r)2 ( (p2 _ fc2^p2g2('2^ _|_ n(g2))(^2^^2j^2 _|_ B'^{k?')) 

2 A(fc 2 )(A(fc 2 ) + A(p 2 ))S(fc 2 ) 


T2Ti{x) 


p‘^q^{l + Il{q^)){A‘^ {k‘^)k^ + B^{k‘^))'^ 

2 B(fc 2 )(B(p 2 ) - B(fc 2 )) 


X TiTiiy) - 


(p2 - fc2)p2g2(i n(g2))(A2(A:2)A:2 + B2(fc2))2 

A{k^)B{P){A{p‘^) - A(fc 2 ))(p 2 + )^ 2 ) 


T 2 T,(y) 


(p2 — k^)p'^q^{l + n(g2))(^2^^2^^2 _|_ ^2^^2^^2 

B{P) 


T2Ti{y) 


+ 


(p2 — k^)p'^q^(^l + n(g2))(^2('^2^^2 _|_ 52 ('^ 2 ^^ 
B(p 2 ) - B(fc 2 ) 

(p 2 — k'^)p‘^q'^{l + n(g 2 ))(^ 2 ('^ 2^^2 _|_ 


T2Ti{y) 

T2Ti{y) 


(60) 


dFA 

dci 


f k'^sinOdkdO f 

J (27r)2 \ 


2l(fc2)(A(A:2) + ^(p2)) 


p 2 g 2 (l + n(g 2 )) 2 (^ 2 (fc 2))!,2 + B‘^{k^)) 

B{P){B{p^) - B(fc 2 )) 


TiTi{z) 


(p2 — k‘^)p^q^{l + n(g2))2(^2^^2^^2 _|. B^{k^)) 

A{k^){A{p^) — A{k^)){p^ + k^) 

2 (p 2 — k^)p^q'^{l + Il{q'^))^{A'^{k^)k^ + B^{k^)) 


T2Ti{z) 


T2Ti{z) 


(61) 


OFb 

dtti 


-! 


k'^sinO dk dO 
(27r)2 


B(fc2) 


g 2 (l + n(g 2 ))(A 2 (fc 2 ))!c 2 + B 2 (fc 2 )) 

B{k^) 


Ti{x) 


+ 


(p 2 — A: 2 )g 2 (l -(- n(g 2 ))(^ 2 ^^ 2^^2 _|_ B^{k^)) 

2 A(fc 2 )(A(fc 2 ) + A{p^))B{e)e 


g 2 (l + n(g2))(A2(fc2))!;2 + B‘^{k‘^)f 
B{P) 


Tiiy) + 


T2Ti{x) 

2A{k^){A{p^)B{P) - A{k^)B{p‘^)) 


g2(l + n(g2))(A2(A:2)fc2 + ^2(^2)) 


?)(y) + 


(p 2 — k^)q‘^(l + Il{q‘^)){A'^ {k^)k^ + B^{k^))^ 

__ 

(p 2 — A: 2 )q 2 (l -|- n(g 2 ))(^ 2 ('^ 2^^2 _|_ B'^{kP‘)) 


T2Ti{y) 

T2Ti{y)Y 


(62) 


dPB , . f khinOdkde I _ ^ _ 

dbi ’ J (27r)2 \(p 2 _A: 2 )(l + n(g 2 ))g 2 (^ 2 (fc 2))^2 + 52 (fc 2 )) 2 ) 

2{A{k^) + A{k^))B{k^) 2B{k^){A{p^)B{k^) - A{k^)B{p^)) 

g2(l + n(g2))(A2(fc2)p + _B2(fc2))2 ^2 + ^^3 _ fc2)g2(l + B{q^)){A^ {k‘^)k‘^ + B^k^))^ ^ 

_ A(fc2) + A(p2) _A(p2)_ 1 

g2(l+n(g2))(A2(fc2)fc2 + 52 (fc 2 )) 2 (p2 _ )^2)g2(l + n(g2))(A2(fc2)fc2 + 52 ()^ 2 )) 2 ^{yI j , 
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(63) 



dFe 

dci 


dFn 

dtti 


f k'^sinOdkdO f 

J (27r)2 


iA{P) + A{p‘^))B{P) 


+ 


(27r)2 1^ g2(l n(g2))2(^2(fc2)p + 52 (^ 2 )) 

{A{p^)B{e) - A{e)B{p^)) 

(p2 - k‘^)q2(l + n(g2))2(^2()!,2)fc2 + ^2(^2)) 


Tiiz) 

T2Ti{z)Y 


f k‘^sinOdkdO j 

J {2ny { 


2A{e)A{p^){A{e)e + A{p^))p^ 
q^{A^{p‘^)p^ + B{p‘^)Y{A^{k‘^)k‘^ + B‘^{k‘^)) 
A{e)A{p^) 


TsTiix) 


TsTiix) 

TsTiix) 


+ 


q‘^{A‘^{p‘^)p‘^ + B{p‘^)){A^(kP‘)kP‘ + B‘^{kP‘)) 

_ A{k^){A{k^) + A{p‘^)) _ 

q^ {A^ {p‘^)p‘^ + B{p‘^)){A^{k‘^)k‘^ + B‘^{k‘^)) 

AA{p^){A{pP‘) — A{k‘^))B{k‘^)B{p^)p^ 

{p‘^ — kP‘)q^ {A^ {p‘^)p‘^ + B{p‘^)Y{A^{kP‘)kP‘ + B‘^(kP‘)) 
2B{e)B{p^) 


T^Tiix) 


+ 


— k‘^)q^{A^{p‘^)p‘^ + B{p‘^)){A^{k‘^)k‘^ + B‘^{k‘^)) 
2A{k‘^)A^{p^){A{pP‘) — A{}A‘))p^ 


{p‘^ — kP‘)q^ {A^ {p‘^)p‘^ + B{p‘^)Y{A^{kP‘)kP‘ + B‘^{kP‘)) 
A{k‘^)A{jP‘) 


TATi{x) 

nTi{x) 


+ 
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where 


_ ;n(p 2 ) - ln(A^)+lnie^) ^ _ ln{K-^)+ln{e^) ^ ^ ^^(^2) _ ln{K^)+ln(e^) 

Zn(A 2 )-in(e 2 ) > 2 /— ln { A ^)- ln { e ^) 

2 2 2 

To avoid the singular points, in integration is split as 

rA np pA 

/ = / + / > (69) 

^ S J £ J p 

where p is the external momentum. In the numerical solution, n-point Gaussian quadrature rule over the internal 
momentum k is used. The number of Gaussian points on two ranges [e^p] and [p. A] are set to Ni = 50 and N 2 = 50. 
In addition, the number of integration over the angle 0 is set to Nq = 40 


in 



3.3 Global convergence 

The newton’s method is sensitive to the initial values of iteration, in this article a method with global convergence 
is presented. The iteration of Newton’s method involves the iteration [18] 


— FiXu)-) 


(70) 


where 



FA{ai,bi,Ci) 

F{x) = I FB{ai,bi,Ci) 
Fu{cii, bi, Ci) 


J{x) = 


dFijx) 

dx^ 


^n+l — ^n+l5 


(71) 


with n indicating the number of iteration. The global converging method modifies the iteration = Xn — A^+i 
by setting [23] 


^n+l — (^^) 

where A G (0,1]. This modification gives us the flexibility in adjusting the step length of iteration. 

The An and F{xn) should approach to 0 in the limit n ^ oo, when x is assumed to converge to the exact value. 
Since F{xn) should approach 0 after sufficient amount of steps of iteration, we expect that our method should 
reduce the absolute value ||F(xn )||2 after every iteration.This can be achieved by carefully adjusting the step length 
of iteration through choosing appropriate A. We set the condition for convergence to be ||F(xn )||2 < where is 
the predetermined absolute tolerance. We also need a relative tolerance as the condition for convergence. Thus 
the final F{xn) should satisfy 

||i^(^n)||2 < T^||i^o(^n)||2 + (73) 

where Fo{xn) is the exact value. 

To avoid possible stagnancy in the iteration, we need two parameter ai and (72, with 0 < ai < (72 < 1, then we 
have the following safeguard against possible stagnancy in the iteration. A in the current trial step should satisfy [23] 


A = 


(7i Ai, A < (7i Ai 

A, (7iAi <c a <c ( 72A1 , 
(72A1, cr2Ai < A 


where Ai is the value in the previous trial step. 


(74) 


3.4 Three-point parabolic model 

In the algorithm, we first set A = 1 to iterate. In a trial step, if the A = 1 fail to reduce the absolute value 
||A„|| 2 , then, we set A = ctiA. When we have two successive trials that fail to reduce the value of ||A„|| 2 , the next 
A can be determined by the three-point parabolic model 


(Ai(F2 - F0)/A2 - A2(F1 - F0)/Ai)/(Ai - A2) 
2((F1 - F0)/Ai - (F2 - F0)/A2 )/(Ai - A 2) 


where 


F0 = F{xn), FI = F{xn - XiAn+i), F2 = F(a;„ - A2A„+i). 


(76) 


4 Numerical results 

For comparison, the unknown functions A{p^), B{p^),Il{p^) with quenched approximation, one-loop contribu¬ 
tion, and the minimal BC vertex (namely the first term in the BC vertex, which introduces the wavefunction 
renormalization dependence in the vertex when compared to bare approximation) are also calculated. 

4.1 Quenched approximation 

In the quenched approximation, we have 

n(g2) = 0, (77) 

i.e. we neglect the fermion loop contribution to the vacuum polarization and adopt BC vertex in the calculation. 
Under this approximation, we calculate the A(p^), 5(p^), n(p^) using the methods stated above. Then we have 






Figure 6: The behavior of A(p^) and B{jA) with Nf — 1 and = 0. 


4.2 One-loop approximation 

In the one-loop, we only consider the one-loop contribution of fermion to the vacuum polarization 

= - / 77^Tr[- -- 7 ^ . (78) 

J [ziry ’ k ' p 

Also we can factor out the tensor structure of the vacuum polarization 

- qfiqu)'n{q^). (79) 

Then we get the contribution of massless fermion to the vacuum polarization scalar 

Af. 

n(g 2 ) = (80) 

The result is shown in Figure 7 




Figure 7: The behavior of AijA) and BijA) with = 1 in the one-loop approximation. 



Figure 8: One-loop contribution to photon polarization. 


4.3 Minimal BC vertex 

In the approximation of the minimal BC vertex, we only assume the first term in the BC vertex. This strategy 
introduces the fermion wavefucntion renormalization dependence in the vertex, while avoid the complication of 
singular structure in BC vertex. 


-N =1 -N =1 




-N=1 



Figure 9: The behavior of and n(p^) in the minimal BC vertex. 


4.4 The complete BC vertex 
4.4.1 The chirally broken phase 

The improved Newton’s method is applied to the calculation of A{jp^)^ and n(g^). The iteration is extremely 

sensitive to initial values. The initial valued can be set to the results of the minimal BC vertex. The Nambu solution 
is shown in Figure 10. 

















-N^= 1 



Figure 10: The behavior of A{p^)^ B{jA) and n(^^) in the BC vertex. 

The above calculations show that the behaviors of dressing functions depends critically on truncation schemes. 
Comparison between the quenched case and the one-loop case shows that the contribution of one-loop vacuum 
polarization changes the behavior of A{p^) significantly, while having little bearings on the behavior of 
After introducing wave-function renormalization in the vertex, although the behavior of A{p‘^) is modified slightly, 
B{p^) remains almost the same. The vacuum polarization 11 (p^) has similar behavior to that of B{p^). When the 
complete BC vertex is considered, again the shape of A{p^) is changed drastically, while B{p^) and 11 (p^) show little 
change. We thus reach the conclusion that the behaviors of B{p^) and 11 (p^)-which represent contributions to the 
fermion and photon self-energy, respectively-are insensitive to truncation schemes, while A(p^)-which contributes 
to the wave-function renormalization-demonstrates close connection to the form of the fermion-photon vertex. We 
also find a model-independent result that although A(p^) shows different behaviors in different approximations, it 
approaches its asymptotic form at p^ 0 and for B{jp^) and 11 (p^) this happens at roughly p^ 10“^. 

4.4.2 The symmetric phase 

The chirally symmetric phase (often referred to as Wigner phase, and the corresponding solution of the gap 
equation is called the Wigner solution) of the DSEs is necessary in the discussion of phase transition problems in 
QEDs and QCD [25h3Q] . This solution can be obtained by setting B(p^) = 0 and it has been found by the authors 
of Ref. [7] in the bare and BCl vertices. However, no solution in the complete BC vertex has been reported before. 
We therefore make an attempt to calculate the Wigner solution in BC vertex adopting the method stated above. 

The power-law form A(p^) = cp^'^ is generally assume in the infrared region. In the case of bare and minimal BC 
vertex (or BCl), self-consistent solutions satisfying DSEs have been found and are related to flavor number by [TT] 

Kbare - ^ ^ + O(^), ( 81 ) 


Id 















and for BCl vertex, two solutions have been found 


I^BCl = 


0.115 


0.044 1 , 


N] 


( 82 ) 


i^BCi = 0.5 — 


0.050 


0.006 

-Nf 


0.028 

-Nf- 


(83) 


However, the complete BC vertex doesn’t admit self-consistent power-law form of solutions. Thus we take the 
empirical extrapolation with k. = 0 . 12 , which can be continuously matched to the numerical solutions with the 
condition 


A{p‘^) A{e^) 


(84) 


where 5 is the infrared cutoff. 

In the ultraviolet region, the solutions restore to the free form. Then we have A{p) = (xl/q. The result 

is shown in Figure 11. 




Figure 11: The Wigner solution of A{p‘^) and n(p^) with = 2 in the BC vertex. 

Due to the complexity of the problem, this solution was not reported in any known literature, and by comparison 
with the results presented in [7], we find that when adopting the complete BC vertex, A{p^) approaches 0 more 
rapidly as p —> 0 , while n(p^) remains insensitive to truncation schemes. 

5 Conclusions 

An improvement is made in the calculation of dressing functions in QED 3 . This method has shown better con¬ 
vergence than other existing methods, especially when applied to the calculation of vertices with singular structure. 
This method can be adopted in the calculation of dressing function in more flavors, and thus determine the flavor 
dependence of quantities such and fermion-antifermion condensate. 

Furthermore, attempts have been made to search for Wigner solution of DSEs in BC vertex. Similarly, we can 
further explore the flavor dependence of Nambu and Wigner solution of DSEs on the basis of this method, the 
results can be applied to problems such as determining the critical point of phase transition. 
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Appendix A. Trace Techniques 

In three dimensional Euclidean space 
tr(/) = 4, 
















tT:{YY) = 4 ( 5 ^'^, 
tr(;5i;^2) = 4pi ■p2, 

tT:{i>ihhi>i) = 4[(Pi ■P2){Pz -Pi) - {Pi ■Pz){P2 -Pi) + {Pi ■P4){P2 -Ps)], 
tr(^i^ 2 i = 0, when n is odd number, 

where 7 ^ is 4-dimensional matrix satisfying 

{7^5 7z^} ~ 2 
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